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Abstract 

A  computer  method  for  objective  computation  of  adiabatic  three- 
dimensional  trajectories  in  the  free  atmosphere  is  presented.     This 
technique  utilizes  a  digital  adaptation  of  the  graphical  iterative  method 
introduced  by  Danielsen  (1961). 

The  program  is  applied  to  a  test  case  and  is  found  to  give  generally- 
satisfactory  results.  Some  difficulties  are  encountered  in  cases  of  very 
low  wind  speeds  and  also  in  regions  of  strong  anticyclonic  wind  shears. 

A  method  for  inclusion  of  diabatic  effects  is  outlined  in  detail. 
This  approach  may  allow  investigation  of  more  complicated  flow  problems 
than  have  been  previously  possible. 


1  .      Introduction 

The  determination  of  three  dimensional  trajectories  by  air  parcels 
in  the  free  atmosphere  has  been  of  interest  to  meteorologists  for  many 
years.     Knowledge  of  parcel  trajectories  is  valuable  in  analysis  of  trace 
substance  transports,  in  the  diagnosis  of  large  scale  systems,  and  also 
can  be  of  importance  in  weather  forecasting. 

Although  the  large  scale  motion  of  the  atmosphere  is  overwhelmingly 
horizontal,  it  has  long  been  recognized  that  neglection  of  the  vertical 
velocity  component  can  produce  appreciable  trajectory  errors  in  as  little 
as  12  hours   (e.g.,   see  Petterssen,   1956).     However,   since  the  large 
scale  vertical  motion  must  be  inferred  rather  than  measured  directly, 
progress  has  been  very  slow  in  developing  satisfactory  trajectory  techni- 
ques . 

Because  diabatic  effects  in  non- precipitating  regions  are  generally 
accepted  to  be  quite  small  in  the  free  atmosphere,  the  concept  of  evalua- 
ting trajectories  on  isentropic  surfaces  has  long  been  an  appealing  one 
(Panofsky,    1947;  Saucier ,    1955;  Petterssen ,   1956).     Many  different 
techniques  have  been  used  previously  --  the  most  straightforward  being 
a  direct  graphical  calculation  from  isotach  and  streamline  analyses 
(Palmen  and  Newton,   1951;  Danielsen,    1959;  Reiter  and  Danielsen, 
1960;  Staley,   1960,    1962;  Reiter,    1963). 

In  any  trajectory  scheme  the  accuracy  of  the  computed  displacement 
tends  to  degenerate  rather  seriously  with  time.     This  is  due  to  the 
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uncertainty  inherent  in  basic  data  which  is  only  available  in  discrete 
form  in  both  space  and  time.     Calculations  of  trajectories  by  the  above 
method  directly  from  streamline  and  isotach  analyses  are  particularly 
difficult  due  to  the  extreme  sensitivity  of  the  calculation  to  errors  in 
the  strsamline  analyses.    For  example,  a  wind  direction  error  of  just 
10     can  produce  a  normal  displacement  from  the  true  trajectory  of  about 
2     latitude  in  just  12  hours  under  average  wind  conditions.    The  effects 
of  errors  of  this  type  and  of  radiosonde  spacing  errors  were  pointed  out 
by  Djuric  (1961). 

In  an  attempt  to  minimize  the  uncertainty  due  to  streamline  errors, 
Danielsen  (1961)  developed  a  more  objective  technique  based  upon  a 
path  integration  of  the  total  energy  equation  expressed  in  isentropic 
coordinates.    This  basic  graphical  technique  has  subsequently  been  used 
with  good  results  by  Danielsen  (1964),  Danielsen,  Bergman,  and 
Paulson  (1962),  Houghton  (1962),  Mahlman  (1965),  and  by  Reiter  and 
Mahlman  (1965a,  b,  c).    This  approach  has  been  used  in  tracing  move- 
ment of  radioactive  debris,  with  computed  results  in  agreement  with 
observed  transports. 
2.    Methods  of  computing  three-dimensional  trajectories 

As  noted  in  the  previous  section  the  energy  method  developed  by 
Danielsen  (1961)  has  proven  to  be  an  effective  technique  for  obtaining 
three  dimensional  trajectories  in  the  free  atmosphere.    A  basic  diffi- 
culty in  employing  this  system  is  the  time  of  computation  involved. 

Initially,  isentropic  surfaces  must  be  generated  containing  analyzed 
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wind  and  Montgomery  stream  function  (c  T  +  gZ)  information.    An  excel- 

P 

lent  treatment  of  the  problems  involved  in  obtaining  Montgomery  stream 
functions  on  isentropic  surfaces  is  given  by  Danielsen  (1959).     Next,  in 
the  evaluation  of  the  trajectories,  a  graphical  iterative  scheme  is 
employed.     This  becomes  very  time  consuming  and  cumbersome  when  a 
large  number  of  trajectories  must  be  determined. 

One  way  to  avoid  this  difficulty  without  loss  of  computational 
accuracy  is  to  perform  the  calculations  with  a  digital  computer  rather 
than  by  graphical  methods.     This,  of  course,  is  advantageous  in  that 
any  such  technique  is  completely  objective,  once  it  is  successfully 
programmed.     Several  previous  investigators  have  developed  computer 
techniques  for  evaluating  three  dimensional  trajectories. 

Nagle  and  Clark  (1966)  employed  a  kinematic  scheme  on  isentropic 
surfaces  using  the  sum     of  the  geostrophic  wind  and  a  divergent  wind 
inferred  from  the  continuity  equation  expressed  in  isentropic  coordinates  . 
This  procedure  may  not  yield  the  most  satisfactory  results  because  Lhe 
divergent  part  of  the  wind  field  inferred  by  this  method  may  be  far  from 
reliable.    Also,  Krishnamurti  (1966a)  showed  that  appreciable  ageostrophic 
wind  components  are  mostly  non-divergent.     This  difficulty  could  possibly 
be  reduced  by  using  the  balance  equation  to  deduce  the  non-divergent 
stream  function. 

Paegle  (1966)  computes  three-dimensional  trajectories  kinematically 
in  a  pressure  coordinate  system  using  a  vertical  motion  (a:)  in  the 


vertical  and  the  sum  of  a  non-divergent  and  divergent  wind  in  the  hori- 
zontal.   The  non-divergent  stream  function  (0)    is  obtained  from  the  non- 
linear form  of  the  balance  equation.    The  velocity  potential  (\)  is  ob- 
tained from  the  vertical  motion  field.    This  vertical  motion  field  is  derived 
by  a  succession  of  iterations  of  the  complete  non-linear  omega  equation 
(Krishnamurti ,   1966b).    This  technique  contains  the  possibility  of  giving 
excellent  results  provided  that  the  derived  three  dimensional  winds  are 
not  overly  sensitive  to  errors  in  the  initial  height  field  and  further,  are 
not  excessively  smoothed  during  the  complex  numerical  procedures.    In 
fact,  part  of  the  objective  of  this  study  is  to  develop  an  alternative 
trajectory  method  to  this  one  in  order  to  provide  more  rigorous  tests 
for  each  approach.    This  is  particularly  important,  since  both  schemes 
are  very  complex  numerically,  and  require  considerable  preparation  of 
input  data . 

In  the  usual  kinematic  trajectory  scheme  one  integrates  the  equations 

d  0  =   J    dt  (1) 

a. 

dX  = — -  dt 

a  cos   0 

where    0    is  the  latitude,    X    the  longitude,    u    and    v    the  horizontal 

scaler  velocity  components,  and    a    is  the  radius  of  the  earth.    Because  the 

u    and    v    velocity  components  are  either  explicitly  or  implicitly  derived 

from  the  analyses  of  streamlines  and  isotachs  ,  errors  in  streamlines 

can  produce  appreciable  errors  in  the  final  trajectory  positions.    In  the 

free  atmosphere  the  wind  shears  normal  to  the  flow  are  generally  much 
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larger  than  those  along  the  flow.  As  a  result  of  this,  streamline  errors 
producing  computed  trajectory  positions  normal  to  the  true  position  are 
especially  serious. 

A  way  to  reduce  the  problems  associated  with  this  difficulty  would 
be  to  devise  an  objective  method  which  "constrains"  the  parcel  from 
deviating  so  far  laterally  from  its  true  position.    The  graphical  method 
proposed  by  Danielsen  (1961)  is  specifically  designed  to  provide  this 
constraint  by  using  the  additional  information  given  by  the  distribution 
of  thermodynamic  variables. 

Because  this  graphical  method  has  produced  highly  reliable  results, 
a  proper  computer  simulation  of  this  technique  should  also  prove  to  be 
reliable.    The  rest  of  this  section  will  discuss  the  concept  behind  this 
method  and  demonstrate  how  it  is  done  on  a  computer  rather  than  graphi- 
cally. 

Since  this  method  is  so  fundamental  to  the  subsequent  development 
of  the  machine  method,  the  basic  concepts  of  this  system  will  be  pre- 
sented here . 

The  horizontal  equation  of  motion  in  isentropic  coordinates  is 

dV  _      _ 

'   =  -  f  kx  V_    -  V^M  +  F  (2) 


dt  ~        2         e 

— » 

where    V      is  the  horizontal  wind  vector,    f    the  Coriolis  parameter,  k 
a  unit  vector  in  the  vertical,    v~   the  two  dimensional  gradient  operator 
on  a  surface  of  constant  potential  temperature    6,    M    the  Montgomery 


stream  function,  and    F    is  the  friction  force. 

Taking  the  dot  product  of  the  horizontal  wind  vector  into  Eq .   (2) 
gives 

df"  ("Vv  =    -  V2    •  V^  +  v2   -  F,  (3) 


the  mechanical  energy  equation  expressed  in    6  coordinates.     Performing  an 

dM      dM 
Eulerian  expansion  on  the  first  term  and  noting  that  -V-  •  V  M  ~  -  — — +  r—    , 

2      e  dte    ate 

V    •  V                             V    •  V 
d  /     2  __2_-\        d6     3   /     2 2\  dM        BM        -        -j 

9  9  6 

The  substantial  derivatives  with    0   subscripts  now  represent  the 
substantial  derivative  computed  without  the  vertical  advection  term. 
If  Eq .   (4)  is  now  integrated  over  the  time  interval  At,    one  obtains 
(1)  (2)  (3) 


2        2 
V^-  V. 
f         l 


tf   SM 


e      t.      e 

(4)  (5) 

•i  ii 

where  the    i    and    f    subscripts  represent  the  initial  and  final  values,  re- 

2       -*        -»  2         2 

spectively,and    V    =  V     '  V    =  u     +  v    .    The  left  hand  side  of  Eq .   (5) 

represents  the  change  of  the  Montgomery  stream  function  over  the  time 
interval  of  the  projection  of  the  true  trajectory  onto  the    8    surface. 
The  first  term  on  the  right  hand  side  of  Eq.   (5)  (term  2)  is  the  change 


of  kinetic  energy  of  the  projected  trajectory.     Term  (3)  represents  the 
integrated  local  M  tendency  along  the  trajectory  path.    Term  (4)  is  the 
change  of  kinetic  energy  on  the  original    0    surface  due  to  non-adiabatic 
motions.     Finally,  term  (5)  is  the  work  done  on  the  mass  of  air  by  friction 

Eq  .   (5)  is  now  in  a  form  such  that  it  can  be  used  to  produce  a  better 
approximation  to  the  final  trajectory  point.    In  general,  this  equation 
must  be  solved  iteratively  since  the  final  value  of    M    determined  from 
it  is  a  function  of  the  true  trajectory  path.     Eq .   (5)  thus  serves  to  provide 
a  constraint  on  the  deviations  of  a  computed  trajectory  normal  to  its  true 
position  resulting  from  streamline  errors.     Since  the  Montgomery  stream 
function  on  a     0   surface  is  exactly  analogous  to  the  geopotential  height 
on  a  pressure  surface,  a  restraint  on  the  permissible  values  of    M 
provides  a  satisfactory  method  for  reducing  the  error  of  a  computed  tra- 
jectory. 

In    graphically  computed  isentropic  trajectories  the  usual  procedure 
is  to  determine  the  first  approximation  to  the  true  trajectory  by  evalua- 
ting Eq  .    (1)  in  the  integrated  form 

S  =  V  At,  (6) 

where    S    is  the  displacement  and    V    is  the  average  wind  speed  along 
the  trajectory  over  the  interval    zit  .    The  information  from  this  initial 
approximation  is  then  used  to  compute  a  value  of    M      from  Eq  .    (5) 
with  terms  4  and  5  omitted.    If  the  computed  value  of    M      from  Eq.    (5) 
differs  from  the  value  of    M     measured    from  the  initial  trajectory,  the 


initial  trajectory  is  then  displaced  laterally  until  the  final  point  on  the 
chart  satisfies  the  computed  value  of    M   .    At  this  point  Eq .   (6)  is  re- 
evaluated along  the  new  path  to  correct  for  possible  changes  in  wind 
speeds  along  this  new  path.    This  new  trajectory  is  then  used  to  re- 
compute the  terms  in  Eq.   (5).    If  the  computed  and  actual  values  agree 
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within  predetermined  limits  (usually   ~  .1  x  10     erg/g  ),  the  trajectory 

is  completed.    If  not,  the  process  is  repeated  until  the  desired  accuracy 
is  reached. 

Even  though,  in  general,  air  movement  in  the  free  atmosphere  is 
neither  frictionless  or  adiabatic,  it  is  commonly  accepted  that  these 
assumptions  are  approximately  valid  over  short  periods  of  time.    According 
to  Kung  (1967)  a  reasonable  value  of  energy  dissipation  in  the  free  atmos- 
phere is  2  .3  erg  g       sec         or     .01x10     erg  g        12  hr  If  this  value 
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is  compared  with  the  typical  values  of  .5  x  10     erg  g        12hr       for  terms 

1,2,  and  3  in  Eq .    (5),  its  neglection  is  clearly  justified.     However,  in 

the  surface  layer  its  effect  becomes  more  important.    The  effect  of 

neglecting  diabatic  motions  is  at  times  more  serious.    Inclusion  of  these 

effects  will  be  discussed  in  detail  in  section  4. 

3.    Machine  adaptation  of  the  energy  integral  method 

The  previous  section  has  shown  how  the  energy  equation  expressed 

in  0    coordinates  can  be  used  to  obtain  more  accurate  three-dimensional 

trajectories  in  the  atmosphere.     Because  of  its  utility  this  system  has 

been  adapted  for  solution  on  the  computer  by  digital  rather  than  graphical 
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methods.     Due  to  the  requirement  of  an  iterative  solution  to  this  problem, 
there  exists  an  almost  unlimited  number  of  approaches  which  could  prove 
successful.     The  method  employed  here  will  be  described  in  some  detail 
to  illustrate  one  solution  to  the  problem  . 

Although  the  longer  range  goal  is  to  apply  the  technique  to  the 
generalized  problem  including  diabatic  flows  (see  section  4)  the  initial 
programming  is  restricted  to  the  case  of  isentropic  flow.     This  produces 
an  appreciable  simplification  in  the  original  programming  logic. 

All  calculations  are  performed  on  a  Cartesian  grid  of  arbitrary  size 
and  spacing.     The  spacing  between  grid  points  is  fixed  in  degrees  lati- 
tude and  longitude.     Because  of  this,  the  actual  spacing  in  the  east- 
west  direction  is  a  function  of  latitude.     Basic  constants  read  in  are 
the  latitude-longitude  spacing  between  grid  points,  the  southern  lati- 
tude and  western  longitude  boundaries,  number  of  grid  points  in  the  x 
and  y  directions,  number  of  levels  in  the  vertical,  and  number  of  time 
steps  to  be  calculated.    The  initial  variables  read  in  are  wind  direction, 
wind  speed,   pressure,  and  Montgomery  stream  function  at  each  point  in 
the  4-dimensional  (x,  y,   0,  t)  grid. 

The  first  computational  step  is  to  convert  the  velocity  field  to  u 

and  v  components.    Also,  the  variation  of  the  grid  point  spacing  along 

the  x-axis  and  of  the  Coriolis  parameter  with  increasing  latitude  is 

determined.     Next,  the  u  and  v  fields  are  linearly  interpolated  to  give 

values  at  times  between  the  12  hour  observation  times  used  for  the 

input  data.    A  simple  difference  is  taken  to  obtain  the  local  12  hour 
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Montgomery  stream  function  tendency.    At  present  this  12  hour  difference 
is  assumed  to  be  representative  of  the  tendency  over  the  entire  12  hours. 
In  future  programs  higher  order  interpolation  schemes  will  be  employed. 
This  will  be  at  considerable  difficulty  due  to  the  additional  data  tabu- 
lations necessary  and  increased  machine  storage  requirements. 

The  first  approximation  to  the  final  trajectory  is  obtained  by  a 
rather  straightforward  kinematic  approach.     Beginning  at  a  grid  point 
the  parcel  is  displaced  for  one  hour  using  the  initial  values  of    u    and 
v    to  determine  the  length  and  direction  of  the  displacement.    At  this 
point  the  time  interpolated    u    and    v    values  are  spatially  interpolated 
from  the  nearest  four  grid  points  to  the  parcel  to  obtain    u  and    v    wind 
components  for  the  parcel  at  its  new  position.    These  new    u    and    v 
values  are  then  used  to  extrapolate  the  parcel  for  another  hour,  at  which 
time  new    u    and    v    components  are  again  found.    This  process  is 
repeated  until  the  12  hour  time  is  reached. 

At  each  point  along  the  trajectory  a  test  is  made  to  determine 
whether  the  parcel  has  left  the  original  grid.    If  the  parcel  is  found 
to  be  outside  this  data  field,  the  trajectory  is  terminated  and  a  new  tra- 
jectory is  begun  at  the  next  starting  grid  point.    If  the  parcel  remains 
within  the  grid,  preparations  are  begun  for  the  iterative  phase  of  the 
computation  using  the  adiabatic,  frictionless  form  of  Eq .   (5).    Values  of 
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kinetic  energy  and  Montgomery  stream  function  are  calculated  at  the 
initial  and  final  points  of  the  12  hour  trajectory.    Also,  the  finite  dif- 
ference form  of  the  local  Montgomery  stream  function  tendency  is  in- 
terpolated at  each  point  and  then  summed  to  get  a  mean  tendency  along 
the  12  hour  trajectory. 

It  might  be  noted  that  this  initial  approximation  to  the  trajectory 
can  provide  a  valid  result  as  long  as  the  original  wind  fields  have  been 
adequately  specified  and  the  time  interpolations  are  satisfactory.    As 
a  separate  problem  it  is  of  interest  to  note  the  differences  between 
these  original,  kinematically  determined  trajectories  and  the  final  ones 
which  also  satisfy  the  energy  equation. 

To  begin  the  iteration  the  difference  is  obtained  between  the  value 
of  M    computed  from  Eq  .    (5)  with  terms   (4)  and  (5)  omitted  and    M 

obtained  from  the  actual  trajectory  (6M) .    If  the  absolute  value  of  this 

7  -1 

difference  is  less  than  .05  x  10     erg  g      ,  the  initial  trajectory  is 

assumed  to  be  correct  and  the  calculation  proceeds  to  a  new  trajectory. 

If  this  is  not  the  case,  the  trajectory  is  then  corrected  in  the  direction  of 

reducing  this  difference.     This  is  accomplished  by  noting  that  in    6 

coordinates , 

9  f    dn  . 

0 

where    V      is  the  total  geostrophic  wind  and    n    is  the  coordinate 

g 

directed  normal  to  the  M  lines.    A  "correction  distance"  can  be  deter- 
mined from  (7)  in  the  finite  difference  form 
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6n  =    -—  (8) 

s 

Note  that    V    ,  the  actual  wind,  is  used  rather  than  the  geostrophic 
s 

wind  here.    This  value  of    6n    obtained  from  Eq  .   (8)  must  then  be  trans- 
lated into    X    and    0  displacements  so  that  the  exact  direction  of  the 
iteration  can  be  determined.     Since  the  correction  must  be  made  normal 
to  the  trajectory  path,  this  is  accomplished  by  the  expressions 

6X=-^-    2  60  =  -6lL    fi   ,  (9) 

a  cos  0    V  a       V 

s  s 

where    u,  v,    and    V      are  evaluated  at  the  final  point  of  the  first  approxi- 
mation to  the  trajectory. 

One  may  note  that  this  technique  does  not  completely  eliminate  the 
difference  between  the  two  values  of    M     unless  the  flow  is  geostrophic. 
However,  the  flow  may  deviate  far  from  the  geostrophic  and  the  method 
will  still  move  the  trajectory  in  the  proper  direction  for  convergence. 
If  the  angle  between  streamlines  and  Montgomery  stream  function  con- 
tours is  greater  than  90    ,  the  method  would  fail.     Since  this  condition 
is  extremely  rare  in  the  free  atmosphere,  the  method  proves  to  be 
generally  effective. 

Once  6X  and  60    have  been  determined  for  the  final  point,  each 
point  along  the  trajectory  is  displaced  a  distance  proportional  to    6X 

and    60    depending  upon  its  time  of  occurrence  along  the  12  hour 

1  9 

trajectory  (eg.,    6X1  hr    =72~6*"6*9hr    =T26^* 

At  this  point  a  new  approximation  to  the  trajectory  has  been  completed 

12 


However,  this  new  approximation  to  the  trajectory  may  be  in  a  somewhat 

different  wind  speed  regime,  thus  necessitating  an  additional  correction. 

This  is  accomplished  by  writing 

6s  =S    .     ,-  S  ,  (10) 

wind       geom 

where    6s    is  the  distance  correction  (to  be  made  parallel  to  the  trajectory) 

S         .  is  the  total  displacement  along  the  new  trajectory  obtained  by 
wind 

summing  1  hour  displacements  using  the  total  wind  at  each  point.    S 

geom 

is  the  distance  obtained  by  summing  the  actual  geometric  distances  be- 
tween the  1  hour  trajectory  points.    Thus    6s    provides  a  correction  for  the 
actual  wind  speed  in  the  vicinity  of  the  new  trajectory.     Note  that  this  is 
only  a  correction  for  the  speed  and  does  not  correct  for  direction.    This 
value  of    6s    is  then  converted  to    X    and    0  displacements  according  to 

BJl— *5-r  *  !#--*!  (11) 

a  cos  0    V  a     V 

s  s 

At  this  point  the  trajectory  has  been  corrected  for  the  Montgomery 
stream  function  change  plus  an  additional  correction  for  distance.    The 
difference  between    M     computed  from  the  adiabatic,  frictionless  form  of 
Eq .   (5)  along  the  new  trajectory  and  the  value  of    M     actually  measured 

from  the  trajectory  end  point  is  then  recomputed.    If  the  absolute  value 

7  -1 

of  this  difference  is  now  less  than  .05  x  10     erg  g      ,  the  trajectory  is 

complete.    If  not,  the  iterative  cycle  is  repeated  until  convergence  is 

reached. 

Normally,  the  iterative  correction  for  a  12  hour  trajectory  is  not 

particularly  large  (usually  less  than  2     latitude).    However,  in  certain 

13 


situations  when  appreciable  accelerations  are  present,  or  the  large 

scale  systems  are  moving  rapidly,  the  corrections  can  be  substantial. 

7  _i 

Fig.    1  is  an  analysis  of  Montgomery  stream  function  (10     erg  g      ) 

for  a  trajectory  test  case  at  0  =  300K,  beginning  22  November  1962, 

1200  GMT.     Pressures  at    6  =  300K  are  given  in  mb. 

A  sample  of  some  of  the  adiabatic  trajectories  obtained  from  the 
machine  method  outlined  above  are  given  in  Fig.  2  for  the  12  hour  period 
beginning  on  22  November  1962,  1200  GMT.    The  dashed  lines  are  the 
kinematic  trajectories  used  as  a  first  approximation  to  the  final  tra- 
jectories.   The  solid  lines  represent  the  12  hour  trajectories  obtained  by 
the  energy  integral  method.    A  single  solid  line  emanating  from  a  point 
means  that  the  kinematic  trajectory  also  satisfied  the  energy  integral 
condition.    The  two  types  of  trajectories  tend  to  be  similar  —  the 
iterated  trajectory  normally  characterized  by  a  small  lateral  correction 
to  the  kinematic  one. 

A  marked  exception  to  this  general  characteristic  may  be  seen  in  the 
lower  left  trajectory  in  Fig.  2.    In  this  case  the  iterated  trajectory  is 
displaced  almost  perpendicularly  from  the  kinematic  trajectory.    The 
cause  of  this  difficulty  may  be  seen  in  Fig.   1,  which  shows  that  this  is 
an  area  of  very  weak  M  gradients  .     Because  the  lateral  distance  cor- 
rection is  made  from  Eq.    (8),  a  region  of  very  low  wind  speeds  may  lead 
to  an  unrealistic  correction  distance.    Further,  because  the  distance 
of  correction  is  relatively  large  for  a  given  value  of    6M,  in  such  cases, 
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the  correction  is  particularly  sensitive  to  analysis  errors  in  the    M    fields. 
Thus  in  this  case,  the  "corrected"  trajectory  becomes  less  accurate  than 
the  kinematic  trajectory.     In  Itiis  program,  and  also  in  graphical  calcu- 
lations, it  has  been  found  to  be  advantageous  to  simply  use  the    kinematic 
trajectory  in  cases  where  the  wind  speed  is  less  than  10  m  sec 

Another  source  of  error  in  the  iterative  calculations  arises  when 
there  exists  a  very  rapid  change  of  wind  direction  perpendicular  to  the 
estimated  trajectory  path.     In  such  cases  the  energy  iterated  trajectory 
direction  may  differ  somewhat  from  the  actual  wind  direction  at  the  final 
point.     This  type  of  error  also  usually  occurs  when  the  wind  speed  it- 
self is  quite  low.    As  above,  the  kinematic  trajectory  is  then  held  to  be 
the  more  accurate  one. 

A  final  difficulty  in  the  iterative  scheme  arises  when  the  wind  shear 
is  strongly  anticyclonic  in  the  region  of  the  trajectory.     In  this  case  the 
changes  of  Montgomery  stream  function  and  of  kinetic  energy  are  nega- 
tively correlated  in  the  iterative  correction  of  an  approximate  trajectory 
using  Eq  .    (5).     For  example,  a  correction  to  increase  the  value  of  M 
will  simultaneously  cause  a  marked  reduction  in  the  final  kinetic  energy. 
A  result  the  computation  may  be  left  as  far  from  convergence  as  in  the  previous 
approximation.     Occasionally  the  iteration  may  even  diverge  from  the  true 
solution.     In  such  cases  the  kinematic   trajectory  is  again  used  as  the 
better  approximation  to  the  true  trajectory. 

In  regions  where  the  above  difficulties  are  not  present,  the  iterative 
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procedure  gives  very  acceptable  results  for  adiabatic  trajectory  calcu- 
lations .    The  inclusion  of  diabatic  effects  produces  considerable 
additional  complication  and  will  be  discussed  in  the  subsequent  section. 
4  .     Incorporation  of  diabatic  effects  into  trajectory  determinations 

Although  the  air  flow  in  the  free  atmosphere  tends  to  be  roughly 
adiabatic  in  character,  diabatic  effects  are  always  present.    At  times 
these  diabatic  effects  may  lead  to  a  true  trajectory  which  is  appreciably 
different  than  that  obtained  by  assuming  isentropic  conditions. 

At  first  it  might  seem  to  be  contradictory  to  consider  diabatic  effects 
within  the  framework  of  a  problem  formulated  in  isentropic  coordinates. 
However,  even  for  diabatic  trajectories,  the    8  coordinate  system  retains 
its  advantage  over    Z    or    p    coordinates  for  several  reasons.    First,  the 
adiabatic  component  of  vertical  motion  is  invariably  larger  than  the 
diabatic  part.    This  implies  that  the  relative  distance  a  parcel  will 
"depart"  from  its  original  isentropic   surface  is  appreciably  smaller 
than  the  departure  from  a    Z    or    p    surface.    Second,  the  diabatic  effect 
of  release  of  latent  heat  of  a  parcel  is  strongly  dependent  upon  the 
cooling  of  an  air  parcel  due  to  adiabatic  ascent.    Thus,  an  accurate 
determination  of  the  three  dimensional  motion  of  a  parcel  is  logically 
obtained  by  utilizing  an  initial  knowledge  of  the  adiabatic  part  of  the 
vertical  motion. 

This  fact  will  now  be  used  to  develop  a  method  for  evaluating  the 
diabatic  components  of  the  motion  as  an  extension  of  the  isentropic 
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method  outlined  in  previous  sections.    From  the  definition  of  a  mean 

value , 

—  a          ■  2  dp  ,        p2  dp       Tj_       ^  2  dp    .         r.  2  dp     _              /n  _, 

a;  At  =   i  -f  dt  ^  f  -7      dt  +   I       -f  dt  +          Tf    dt              12 

■J  dt           J  dt    ,            ^     dt,            J      dt.. 

t1                     t1  ad  tl       lr            t:      lh 

where    10=  dp/dt,  dp/dt        is  the  adiabatic  part  of    co,  dp/dt       is  the 

ad  lr 

contribution  to    go  from  long  wave  cooling,  and    dp/dt..     is  the    co  due 

ih 

to  release  of  latent  heat  of  condensation. 

However,  the  second  and  third  terms  on  the  right  hand  side  of 
Eq .    (12)  represent  a  pressure  change  due  only  to  a  change  of  potential 
temperature  of  the  parcel.    Thus  one  may  write, 

*At  =  Apad+f  A  9lr  +  H   Mh  <13» 

where    bp/b  9  is  the  measure  of  the  mean  static  stability  along  the  tra- 
jectory and     A8    is  the  change  of  potential  temperature.    Thus,  in  this 
model  the  total  vertical  displacement  of  an  air  parcel  is  due  to  a  combi- 
nation of  three  effects  ,  only  the  first  of  which  is  available  from  the 
method  described  previously.    The  term  in  Eq .    (13)  involving    A0.      can 
be  extremely  complex  to  evaluate,   since  it,  in  general,  includes  long 
wave  radiative  effects  of  water  vapor,  ozone  and  carbon  dioxide.    How- 
ever, the  contribution  of  ozone  is  very  small  in  the  troposphere.    Also, 
according  to  Manabe  and  Mo  Her  (1961)  the  contribution  of  carbon 
dioxide  is  considerably  smaller  than  that  of  water  vapor. 

Making  use  of  the  above  fact,  Danard  (1968)  has  developed  a 
technique  for  computing  the  long  wave  cooling  due  to  an  arbitrary  cloud 
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and  moisture  distribution.    Danard's  model  has  been  adapted  for  appli- 
cation to  this  problem.     Fig.  3  shows  diabatic  cooling  rates  (deg .  C/12  hr) 
computed  from  Danard's  model  and  interpolated  to  the    6=  300K  surface 
for  22  November  1200  GMT.     This  figure  shows  that  the  spatial  variation 
of  diabatic  cooling  on  a    6    surface  is  surprisingly  large.    Cooling  rates 
in  excess  of  2  deg/12  hr  are  observed  in  the  northeast  corner.    These 
high  values  are  due  to  radiation  from  the  top  of  the  cloud  cover  on  the 
east  side  of  the  trough  observed  in  Fig.   1.     Note  also  the  low  cooling 
rates  in  the  central  part  of  Fig.   3.     This  area  is  within  and  behind  the 
base  of  the  trough  and  has  been  shown  to  be  a  large  mass  of  descending 
dry  stratospheric  air  in  a  previous  paper  (Reiter  and  Mahlman,   1965c)  . 
Another  reason  for  the  pronounced  variability  of  the  cooling  rates  is  due 
to  the  large  slope  of  the    6=  300K  surface.    In  this  case  the  pressure 
varies  from  920  to  400  mb,  thus  leading  to  large  variations  of  height 
and  temperature  over  the  map.    It  may  be  seen,  however,  that  this 
method  contains  the  possibility  of  gaining  a  more  adequate  determination 
of  the  diabatic  cooling  of  a  parcel  due  to  long  wave  radiation.    The  most 
straightforward  way  to  include  this  in  a  trajectory  calculation  would  be 
to  simply  use  the  average  of  the  cooling  rates  computed  at  the  initial 
and  final  points  of  the  adiabatic    trajectory. 

The  last  term  in  Eq  .    (13)  represents  the  contribution  to  the  total 
pressure  change  of  a  parcel  due  to  release  of  latent  heat  of  condensation. 
This,  of  course,  is  a  heating  effect  and  tends  to  balance  the  diabatic 

cooling  due  to  long  wave  radiation.     The  heating  of  an  air  parcel  due 
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to  condensation  of  water  vapor  may  be  written  with  sufficient  approxi- 
mation as 

f  =¥f "?L-ir  6(x'y<9't)  (14) 

lh  lh 

dH 
where   - —  is  the  heating  due  to  release  of  latent  heat,  L  the  latent  heat 

dtih 

of  vaporization,     Sq  /dp    the  rate  of  change  of  saturation  specific 

humidity  with  respect  to  pressure  and    6(x,y,0,t)  is  either  one  or  zero 

depending  upon    whether  the  point  x,  y,  6,  t  is  saturated  or  not. 

To  obtain  the  actual  value  of  heating  of  a  parcel  over  a  12  hour 

time  interval  due  to  latent  heat  release,  Eq  .   (14)  must  be  integrated 

over  At    ,  the  interval  of  time  which  the  parcel  is  saturated.    This 

value  of    At    is  determined  by  first  determining  the  initial  value  of    q 

for  the  parcel  and  then  finding  the  pressure  (p       )  at  the    8    surface 

sat 

for  which  the  value  of    q      equals  the  initial  value  of    q.    From  this  one 

s 


may  define 


Ap=p.-psat  (15) 


where    Ap     is  the  pressure  change  a  parcel  must  undergo  to  experience 
saturation,  and    p.    is  the  initial  pressure.    Eq  .   (15)  may  then  be  used 


to  form  the  expression 


At'  =  tf  -  t'  =  12  (l  - ^ ^         (i6) 


*>w!K> 


where,  if    At     is  less  than  zero,  it  is  set  to  be  zero  for  the  computation. 
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Here,     (  Ap       +  T^,  A  9,  )  represents  the  total  12  hour  pressure  change 
ad      S  6         lr 

experienced  by  the  parcel  due  to  the  effects  of  adiabatic  motion  and  long 
wave  radiation. 

One  may  now  integrate  Eq .   (14)  over  the  time  interval    At     to  obtain 

A9m  =  -^f¥  ?<"'  or  (16a) 

t       p 


y  co  ?  x  Ba 

A  -1000_L     ^  it^i001L    ^  ^   flt 

lh           c  x                          c             x    £p 

P  P                               P         P 


where    x  =  0.2857.    The  approximation  in  Eq  .   (16b)  is  justified  because 

all  terms  involving  deviations  from  the  mean  in  the  expansion  of 

go   9 q  /Sp/p      are  more  than  one  order  of  magnitude  smaller  than  the  term 
s 

retained. 

If  the  expression  for   A6L     in  Eq .   (16b)  is  substituted  into  Eq .   (13), 
one  obtains 

i^Ap.n-  .gJH&.Sali.'i.^     (17) 

1  ad       B  8        lh      §8        c  x  3p 

P  P 

Note  that  the    cc    values  obtained  from  the  adiabatic  trajectory  and  the 
long  wave  cooling  effect  are  used  as  the  initial  approximation  to  the 
true    co    field  required  in  the  third  term.    The  subscript  on   W,    denotes  that 
this  is  the  first  approximation  to  the  total  pressure  change  by  the  parcel. 
Eq  .    (17)  can  now  be  solved  iteratively  for  the  final  mean    co  by  writing 

To-        At  =  AP        -^A0     _|2l00A    ^_^At'  (18) 

n+1  pad         S6        lr      S8        c  x      3p 

P         P 
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where  the    n    subscript  indicates  the  number  of  iterations .    This  iteration 

can  be  continued  until     I  to      .   -  co      I  <  €,  where    e    is  some  acceptable 

n+1  n 

error.    The  mean  vertical  motion  over  the  12  hour  period  has  now  been 
obtained.    Now  that  the  total  vertical  displacement  has  been  obtained, 
the  total  diabatic  heating  of  a  parcel  is  given  by  the  expression 


Ae=Ae    .iooo^L    %+L    ilsAt. 

lr  c  p>-        dp 

where    A8    determined  along  the  trajectory.    An  interpolation  in    8    is 
utilized  to  find  wind  speeds  at  each  new  potential  temperature  along  the 
path.    Utilizing  these  new  values,  a  correction  is  then  made  for  distance 
according  to  Eq  .   (10). 

To  begin  preparation  for  the  lateral  correction  of  the  diabatic  trajectory, 
A0    from  Eq.   (19)  is  substituted  into  term  (4)  of  Eq .   (5)  written  as 


6i 
The  vertical  derivative  of  the  kinetic  energy  is  evaluated  along  the  tra- 
jectory by  a  centered  finite  difference  of  the  kinetic  energies  of  the  levels 
above  and  below  the  parcel  at  a  given  point  along  the  trajectory. 

The  computed  value  from  Eq .   (2  0)  is  then  substituted  into  Eq.   (5)  to 
obtain  a  new  computed  value  of  Mf  for  the  trajectory.    A  lateral  correction 
is  then  made  for  the  trajectory  in  the  same  manner  as  done  previously  for 

the  adiabatic  case.    This  correction  is  normally  small  because  the  numerical 

7  -1 

value  obtained  from  Eq .   (2)  is  nearly  always  less  than  .2  x  10     erg  g 

considerably  smaller  than  the  usual  magnitude  of  terms  (1),   (2),  or  (3)  in 
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Eq.    (5). 

If  this  correction  is  small  enough  that  the  computed  and  observed 
values  of    M,    agree  within  the  predetermined  limits,  the  trajectory  is 
complete  and  simultaneously  satisfies  the  energy  equation  and  the 
diabatic  effects . 

Provided  that  a  correction  is  necessary  in  order  to  satisfy  Eq  .   (5), 
the  next  step  is  to  provide  an  additional  correction  for  the  change  of 
distance  along  the  trajectory  due  to  effects  of  horizontal  wind  shear. 
This  corrected  trajectory  is  then  checked  by  a  recomputation  of  Eq  .   (5). 
If  the  computed  and  observed  values  of    M      (projected  from  the  diabatic 
trajectory  onto  the  original    8    surface)  do  not  agree  within  the  predeter- 
mined limits,  this  process  is  repeated  until  the  convergence  criterion  is 
satisfied . 

However,  even  though  Eqs.   (5)  and  (6)  are  now  satisfied  for  this  new 
trajectory,  the  possibility  remains  that  the  computed  vertical  displacement 
is  again  in  error,   since  a  path  correction  has  been  made.    A  new  check  is 
accomplished  by  recycling  through  the  process  outlined  from  Eq .   (14) 
through  Eq  .   (18).    If    u;  is  sufficiently  close  to  the  previous  value 

of    a*  obtained  from  the  first  cycle,  the  trajectory  is  complete.    If  these 

two  computed  values  differ  appreciably ,  the  cyclic  procedure  beginning  with 
Eq  .   (19)  is  repeated  until  one  of  the  exit  criteria  is  satisfied. 

At  the  present  time  little  quantitative  information  is  available  as  to 
the  number  of  cycles  required  for  the  computed  trajectory  to  satisfy  the 
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convergence  criteria  .    In  a  set  of  graphical  computations  of  diabatic  effects 
on  trajectories  by  the  author  in  a  previous  publication  (Reiter  et.  al.  , 
1965),  only  one  to  three  iterations  were  normally  required.    In  that  com- 
putation, conservation  of  equivalent  potential  temperature  was  imposed 
as  an  additional  criterion  to  be  satisfied  along  the  trajectory. 

Although  there  are  no  fundamental  difficulties  in  programming  the 
inclusion  of  diabatic  effects  for  machine  computation,  in  actual  practice 
this  will  prove  to  be  very  difficult  and  time  consuming.    Further,  it  should 
be  noted  that  this  method  will  be  incapable  of  determining  trajectories  in 
areas  dominated  by  sub-synoptic  scale  motions. 
5 .      Summary 

A  working  computer  method  has  been  outlined  for  objective  computation 
of  isentropic  trajectories  in  the  free  atmosphere  using  a  modification  of 
the  scheme  originally  derived  by  Danielsen  (1961).    This  method  has  been 
found  to  give  highly  acceptable  results  in  most  cases.    The  method  was 
found  to  give  poor  results  in  either  areas  of  very  low  wind  speeds  or 
regions  characterized  by  strong  anticyclonic  wind  shears. 

A  general  technique  for  expanding  this  method  for  inclusion  of  diabatic 

flows  has  been  presented  in  some  detail.    This  should  enable  investigation 

of  more  general  three-dimensional  flow  problems  than  have  been  capable 

of  solution  up  to  the  present  time.    The  required  expansions  of  the 

isentropic  trajectory  program  to  include  diabatic  effects  are  currently  in 

progress . 
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